module scvt_mod

  use random_number_mod
  use container
  use string
  use flogger
  use const_mod
  use sphere_geometry_mod
  use delaunay_voronoi_mod

  implicit none

  private

  public scvt_set_quasiuniform_density
  public scvt_run
  public scvt_final

  real(8), allocatable :: rho(:,:)
  real(8), allocatable :: rho_lon_bnd(:)
  real(8), allocatable :: rho_lat_bnd(:)

contains

  subroutine scvt_set_quasiuniform_density(nlon, nlat)

    integer, intent(in) :: nlon
    integer, intent(in) :: nlat

    real(8) dlon, dlat
    integer i, j

    call scvt_final()

    allocate(rho(nlon,nlat)); rho = 1
    allocate(rho_lon_bnd(nlon+1))
    allocate(rho_lat_bnd(nlat+1))

    dlon = pi2 / nlon
    do i = 1, nlon + 1
      rho_lon_bnd(i) = dlon * (i-1)
    end do
    dlat = pi / nlat
    do j = 1, nlat + 1
      rho_lat_bnd(j) = dlat * (j-1) - pi05
    end do

  end subroutine scvt_set_quasiuniform_density

  subroutine scvt_final()

    if (allocated(rho        )) deallocate(rho        )
    if (allocated(rho_lon_bnd)) deallocate(rho_lon_bnd)
    if (allocated(rho_lat_bnd)) deallocate(rho_lat_bnd)

  end subroutine scvt_final

  subroutine scvt_run(max_iter)

    integer, intent(in) :: max_iter

    integer iter, iDVT
    type(voronoi_vertex_type), pointer :: VVT
    type(delaunay_vertex_type), pointer :: DVT, DVT1, DVT2
    type(linked_list_iterator_type) it_DVT

    real(8), allocatable :: rho0(:)
    real(8), allocatable :: lon0(:)
    real(8), allocatable :: lat0(:)
    real(8) rho0_VC, x0(3), x1(3)
    real(8) area, wgt, dist
    real(8) Linf, L2
    logical converged
    real(8), parameter :: eps = 1.0d-10

    allocate(rho0(global_DVTs%size))
    allocate(lon0 (global_DVTs%size))
    allocate(lat0 (global_DVTs%size))

    converged = .false.
    do iter = 1, max_iter
      call log_notice('Iteration ' // to_str(iter) // ' of ' // to_str(max_iter) // ' ...')
      Linf = 0; L2 = 0
      ! Calculate density on Delaunay vertices.
      do iDVT = 1, global_DVTs%size
        DVT => get_DVT(global_DVTs%value_at(iDVT))
        rho0(DVT%id) = get_density(DVT%lon, DVT%lat)
      end do
      ! Calcuate the centroidal point of each Voronoi cell subject to the given density.
      do iDVT = 1, global_DVTs%size
        DVT => get_DVT(global_DVTs%value_at(iDVT))
        it_DVT = linked_list_iterator(DVT%linkDVTs)
        x0 = 0; wgt = 0
        do while (.not. it_DVT%ended(cyclic=.false.))
          DVT1 => get_DVT(it_DVT%value)
          DVT2 => get_DVT(it_DVT%next_item%value)
          x1(1) = DVT%x * rho0(DVT%id) + DVT1%x * rho0(DVT1%id) + DVT2%x * rho0(DVT2%id)
          x1(2) = DVT%y * rho0(DVT%id) + DVT1%y * rho0(DVT1%id) + DVT2%y * rho0(DVT2%id)
          x1(3) = DVT%z * rho0(DVT%id) + DVT1%z * rho0(DVT1%id) + DVT2%z * rho0(DVT2%id)
          area = calc_area([DVT%x, DVT1%x, DVT2%x], [DVT%y, DVT1%y, DVT2%y], [DVT%z, DVT1%z, DVT2%z])
          x0 = x0 + area * x1
          wgt = wgt + area * (rho0(DVT%id) + rho0(DVT1%id) + rho0(DVT2%id))
          call it_DVT%next()
        end do
        x0 = norm_vector(x0 / wgt)
        call inverse_cartesian_transform(lon0(iDVT), lat0(iDVT), x0(1), x0(2), x0(3))
        dist = calc_distance(lon0(iDVT), lat0(iDVT), DVT%lon, DVT%lat) / radius
        L2 = L2 + dist**2
        if (dist > Linf) Linf = dist
      end do
      L2 = sqrt(L2 / global_DVTs%size)
      print *, L2
      if (L2 < eps) then
        converged = .true.
      else
        do iDVT = 1, global_DVTs%size
          DVT => get_DVT(global_DVTs%value_at(iDVT))
          call DVT%set(lon=lon0(iDVT), lat=lat0(iDVT))
        end do
        call delaunay_voronoi_reinit()
        call delaunay_triangulation(init_three_vertices=[11,5,2], mute=.true.)
      end if
      if (converged) exit
    end do
    call voronoi_diagram()

    deallocate(rho0, lon0, lat0)

  end subroutine scvt_run

  real(8) function get_density(lon, lat) result(res)

    real(8), intent(in) :: lon
    real(8), intent(in) :: lat

    integer i, j
    logical found

    found = .false.
    do i = 1, size(rho_lon_bnd) - 1
      if (lon >= rho_lon_bnd(i) .and. lon <= rho_lon_bnd(i+1)) then
        found = .true.
        exit
      end if
    end do
    if (.not. found) then
      print *, lon * deg
      print *, rho_lon_bnd * deg
      stop 111
    end if

    found = .false.
    do j = 1, size(rho_lat_bnd) - 1
      if (lat >= rho_lat_bnd(j) .and. lat <= rho_lat_bnd(j+1)) then
        found = .true.
        exit
      end if
    end do
    if (.not. found) then
      print *, lat * deg
      print *, rho_lat_bnd * deg
      stop 222
    end if

    res = rho(i,j)

  end function get_density

end module scvt_mod
